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N omenclat ure 


E = elastic modulus (GPa) 
v — Poisson’s ratio (mm/mm) 
a y = yield stress (MPa) 

STRI65 = quadratic triangular shell element in ABAQUS [1] 

S8R = quadratic reduced-integration shell element in ABAQUS [1] 

C3D10 = quadratic tetrahedral elements in ABAQUS [1] 

C3D15 = quadratic wedge element in ABAQUS [1] 

C3D20(R) = quadratic brick element (reduced-integration) in ABAQUS [1] 
a = crack length (cm) 
n = number of cracks in discrete-source 

0 = orientation of discrete-source damage, angle between positive x axis and nearest crack 
d x = distance from middle stiffener to center of discrete-source damage (cm) 

Ki,Kii,Kiii — mode I, II, and III plane strain stress intensity factors (MPay^) 

KicKiic = plane strain fracture toughness for modes I and II (MPay^) 

Pmax = damage-dependent allowable traction, residual strength (MPa) 

P = applied traction (MPa) 

MSE = mean squared error as defined by Eq. (2) 
c v = correlation coefficient as defined by Eq. (3) 

CTD = magnitude of relative displacement between upper and lower fracture surfaces (mm) 
CTD crit = critical value of CTD (mm) 

CTD i , CTD a, CTD hi = opening, in-plane sliding, out-of-plane shearing components of CTD 
(mm) 

d = fixed characteristic distance behind crack front where CTD is monitored (mm) 
da = crack extension (mm) 

(n) = script to denote mesh at n th crack increment 

(n + 1) = script to denote mesh at (n + l) th crack increment 
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I. Introduction 


Resilient aircraft control involves adaptive responses to off-nominal flight conditions, including 
the incurrence of structural discrete-source damage during flight. Discrete- source damage is typically 
manifested as a result of a structural impact event, including hail- and birdstrike. In 2003, an Airbus 
A300 operated by DHL was struck by a surface-to-air missile after takeoff from Baghdad, Iraq, 
causing discrete-source damage to crucial control surfaces of the left wing [2]. In 2008, a Boeing 
747-438 operated by Qantas Airways incurred in-flight structural damage to the fuselage and right 
wing leading edge following the failure of an onboard oxygen cylinder [3]. Although the aircraft 
landed safely in both cases, these examples motivate a need for more resilient, adaptive control 
system responses. 

In these types of cases, problems associated with in-flight discrete-source damage, for example 
inability to sustain original design loads, can be exacerbated by crack propagation from damaged 
regions. To avoid unstable crack propagation, load levels must be maintained below a reduced 
load-carrying capacity, or residual strength , of damaged flight structures. Adaptive control system 
responses might include automatic adjustment of certain flight parameters (e.g. velocity, maximum 
acceleration) to accommodate structural residual strength. This accommodation implies that accu- 
rate residual strength predictions of flight structures with complex damage configurations be made 
in real time , during flight ; this capability currently does not exist for commercial aviation. 

Challenges to developing an adaptive response technology include accurately predicting residual 
strength of discrete-source damaged structures both offline (i.e. during control system design) and 
online (i.e. in real time onboard the aircraft). In the offline context, researchers have developed 
various tools for determining residual strength of thin, damaged metallic structures using elastic- 
plastic fracture mechanics (EPFM)-based numerical methods. For example, two common finite 
element (FE) modeling techniques involve nodal release and adaptive remeshing. Both techniques 
represent cracks geometrically [4]. The former, however, prescribes possible crack trajectories, which 
introduces inherent mesh dependencies into fracture simulations and limits generality of crack path 
predictions. Nodal release techniques have been used in 2D [5-12] and in 3D [13-15] for studying 
crack growth parameters and predicting residual strength of structures where the crack path was 
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known a priori and where mesh refinement along the crack path sufficiently characterized growth 
increments. Adaptive remeshing techniques avoid such mesh dependencies and enable simulation 
of arbitrary crack propagation using evolutionary models or criteria [16-20]. Adaptive remeshing 
techniques have been implemented in both 2D [21] and in 3D [22, 23]. Of the described techniques, 
3D, adaptively remeshed, elastic-plastic tearing simulations provide the most general prediction 
capabilities for crack growth and residual strength. 

It is infeasible to perform a rigorous and computationally intensive crack growth simulation 
within the possible short time span following a discrete-source damage event. Thus, an approxima- 
tion, or surrogate model , is needed for making online predictions of residual strength. Queipo et al. 
provided a complete description of surrogate modeling development and optimization [24]. With 
regard to surrogate construction, they described both parametric (e.g. polynomial regression and 
Kriging) and nonparametric (e.g. radial basis functions) approaches. In nonparametric approaches, 
a global functional form relating system input to system response is not assumed. 

Artificial neural networks (NNs) are a nonparametric surrogate modeling approach and are 
trained to infer a nonlinear mapping from system input to system response, or output. The reader is 
referred to [25] for an extensive methodology overview of the most widely used types of NN. Different 
types of NNs have been applied extensively for damage detection [26-32] and, to a much lesser extent, 
for damage assessment. Ouenes et al. employed a NN methodology to predict fracture indicators 
(e.g. density of fractures) in naturally fractured rock reservoirs as a function of various geological 
and geophysical data [33]. Pidaparti et al. employed a NN to predict residual strength and corrosion 
rate of aging aircraft panels with collinear multi-site damage by training with experimental results 
and validating with both experimental results and analytical solutions [34]. Recently, Mohanty et 
al. used a Gaussian process (GP) approach to predict fatigue crack growth in aluminum 2024-T351 
specimens by training two distinct models, one presented with experimental load parameters as 
input and another presented with piezoelectric sensor signals as input [35]. In that work, Mohanty 
et al. used observed fatigue crack lengths and growth rates as known output for training each model. 

Alternatively, NNs can be trained using results from numerical experiments, or simulations [36]. 
For example, Sankararaman et al. recently used linear-elastic fracture parameters computed from 


4 



FE analyses to train a GP model as part of a method to statistically infer equivalent initial flaw 
size in fatigue applications [37]. High-fidelity numerical simulations can provide training data when 
analytically- and experimentally-derived data are limited due either to a lack of generally applicable 
analytical solutions or to prohibitive costs of obtaining sufficient experimental data. 

The purpose of the work presented here is two-fold: (1) to illustrate a methodology for creating 
a surrogate model as a real-time residual strength prediction tool and (2) to describe and validate 
numerical tools for making accurate residual strength predictions offline using fully 3D, elastic- 
plastic, FE-based crack growth simulations. The high-fidelity, more computationally expensive tools 
described in (2) can provide training data that, when coupled with the surrogate model methodology 
described in (1), can be used in the design of adaptive response technology. 

Consistent with our two-fold purpose, this paper is divided into two primary sections. Section II 
illustrates the methodology for developing a surrogate model (in particular, a NN) that predicts 
residual strength as a function of discrete-source damage parameters. The methodology is illustrated 
using a relatively simple proof-of-concept example. The procedure for gathering training data 
is described in II A and II B. Because an implementation-ready NN is beyond the scope of this 
paper, training data for the proof-of-concept example relies on reduced-order residual strength 
approximations. After collecting training data, a simple NN is constructed in II C by optimizing 
certain performance parameters. Finally, a sensitivity study is conducted in II D to understand the 
effect of each damage parameter on predicted residual strength specifically for the proof-of-concept 
structure. 

Section III improves upon offline residual strength prediction tools used in Section II by simulat- 
ing 3D, elastic-plastic tearing. The tools provide more general crack growth simulation capabilities 
and can be used to generate accurate residual strength training data. A relatively large, integrally- 
stiffened panel (ISP) that exhibits crack branching is simulated in IIIC to validate the tools. 

Results and discussions from the NN proof-of-concept example and from the elastic-plastic tear- 
ing simulation are provided in each respective section. Section IV offers a summary and conclusions 
for the entirety of this work. 
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II. Neural Network Development and Methodology 


This section describes the development of a surrogate model for predicting residual strength of 
discrete-source damaged aircraft structures in real time. A global functional form is not assumed 
for the nonlinear relationship between residual strength and the damage parameters influencing it; 
thus, a nonparametric surrogate model is developed. In particular, a supervised NN is considered 
due to rapid prediction capabilities amenable to real-time applications. In Fig. 1, the upper dashed 
region shows the generalized procedure for developing a NN (surrogate model) that predicts residual 
strength as a function of parameterized discrete-source damage. The lower dashed region shows the 
functionality of the NN (surrogate model) in a real-time context. 



Fig. 1 Upper dashed box illustrates a general approach for developing a surrogate model to 
predict residual strength of damaged structures. Lower dashed box illustrates how the sur- 
rogate model would function onboard an aircraft for predicting residual strength of damaged 
structures in real time. 


The first step in this type of surrogate model development is typically referred to as design of 
experiment (DOE) [24] and involves obtaining data points that will be used to train and test the 
NN. The DOE should be based on the intended application of the NN. For example, if the NN is 
intended to provide residual strength predictions in terms of maximum allowable bending moment 
in a damaged aircraft wing, then the data points should be gathered using an appropriate wing 
structure with applied boundary conditions of interest. Each data point includes sampled input 
variable(s) and corresponding known system response(s), called target output. Once the NN has 
been trained to map given input to target output, it becomes a useful tool for predicting system 
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response when presented with new input that is within the training range but does not necessarily 
correspond to data points used for training. 

To illustrate the methodology, a simple NN is developed using a representative wing structure 
and reduced-order (linear-elastic) approximations for predicting residual strength. The representa- 
tive wing structure is a 61.0 x 91.4 cm 2 integrally-stiffened panel (ISP) with three blade stiffeners 
each 5.1 cm in height, as shown in Fig. 2. The ISP skin and stiffeners are 2.3 mm thick. The panel 
is modeled as linear-elastic with E = 71.0 GPa and v — 0.33, similar to values for a 2XXX series 
lower wing skin aluminum alloy (AA). 


P 



G* I I I ,1 5.1 


(thickness of skin and stiffeners = 0.23 ) 


Fig. 2 Schematic of ISP model with dimensions similar to those used in [38]. Plan view (top) 
and cross-section showing integral blade stiffeners (bottom). A damage-containing region is 
modeled using 3D solid elements (enclosed in shell-solid boundary) while remainder of panel 
is modeled with shell elements. All dimensions in cm. 


Multiple FE models of the uncracked panel are constructed using ABAQUS® [1]. A shell- 
solid modeling technique is employed, where each panel is modeled using 3D solid elements in a 
region that will contain damage and shell elements elsewhere, as depicted schematically in Fig. 2. 
In this way, 3D constraint is inherently captured along crack fronts using fully 3D solid elements, 
while shell elements help maintain a level of computational efficiency yet are able to capture out- 
of-plane deformation and possible buckling. The shell and solid element regions are joined using 
a coupling constraint, whereby resultant forces and moments acting at shell edge nodes on the 
shell-solid boundary are distributed as forces acting at nodes located in a region of influence on the 
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solid surface of the shell-solid boundary. A mesh refinement study is carried out to ensure adequate 
discretization of the panel models. Uncracked panels are modeled using approximately 50 STRI65 , 
2000 S8R, and between 1800 and 17,300 C3D20R elements, depending on the size of the damaged 
region. Boundary conditions for the ISP models are defined to emulate tensile loading conditions 
for a region of the lower wing surface and are shown schematically in Fig. 2. 

A supplementary study was carried out to determine shell-solid boundary effects on nearby 
crack fronts in order to minimize the size of the solid region without affecting stress intensity factors 
(SIFs) computed along nearby crack fronts. Maintaining fracture parameter accuracy is especially 
important since fracture parameters are used to predict structural residual strength (described in 
II B), which is in turn used to train the NN (described in II C). The supplementary study considered 
a 61.0 x 91.4 cm 2 unstiffened panel of the same (linear-elastic) material and thickness as the ISP 
described above. The panel had a single, 12.7 cm long, centrally-located through-crack oriented in 
the x direction (normal to applied tensile load). Both tensile and bending conditions were considered 
in the study. The panel was modeled entirely with shell elements except for a region containing the 
crack, which was modeled with 3D solid elements. All model parameters remained constant while 
varying the size (both in-plane dimensions) of the square-shaped solid region, therefore varying the 
distance from the shell-solid boundary to the crack front. The size of the solid region was initially 
slightly larger than the length of the crack and was increased until computed SIFs converged. Results 
from the supplementary study indicated that for a static, linear-elastic crack, the distance from shell- 
solid boundary to nearest crack front should be no less than 25% of the crack length. This ensures 
that the shell-solid boundary has negligible effect on computed SIFs. The same rule-of-thumb is 
applied to the example ISP models described in the NN study. 

The following sections describe the generally applicable methodology for developing a NN as a 
real-time residual strength prediction tool. 

A. Input Variables: Discrete-source Damage Parameters 

Discrete-source damage in this work is represented by a symmetric, star-shaped, array of equi- 
length cracks, as depicted in Fig. 3(b). This representation of discrete-source damage is motivated 
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by observations of petaling caused by penetration damage to thin metallic structures, see Fig. 3(a). 


If all of the cracks in the star-shaped array of Fig. 3(b) separate under load (i.e. there are no crack 


closure effects), then the cracked region transfers no load and effectively represents a circular hole 


with petaling edges, similar to that shown in Fig. 3(a). The damage representation is parameterized 


by the four variables n, a, d x , and #, which are postulated to influence residual strength of the ISP. 



Reprinted with permission of 
ASM International®. 
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Fig. 3 (a) Petaling on the reverse side of a metallic sheet subject to explosive, discrete-source 
damage [39]. (b) Schematic showing the representation and parameterization of discrete- 

source in the NN example described in this work. 


The sample space of damage configurations is defined by a range of values for each parameter. 
Ranges can be specified based on accident reports, photographic evidence, potential structural 
threats, design specifications, and so forth. Inherently, the NN predictions are valid only for input 
parameter values within the range of training data. Thus, it is necessary to define the sample space 
based on the particular NN application. In the example NN, ranges for each damage parameter 
are limited to some extent by the ISP geometry. Each range is given in Table 1. The parameter 
n, takes integer values ranging from two to six. The range of 0 depends on n due to the definition 
of orientation and the symmetry of the star-shaped configuration. The range of a is defined in 
terms of ISP bay width, from 1/8 * bay width to 1/4 * bay width. Due to symmetry of the ISP 
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model, the parameter d x ranges from 10.2 cm (damage centered in mid-bay) to 0 (damage centered 
at middle stiffener). If the damage is located such that the damage-containing, solid FE region 
overlaps anywhere with the middle stiffener, the stiffener is assumed to be severed in the damaged 
region and is modeled explicitly as such. 


Table 1 Range of values associated with each damage parameter in the example NN. 


Damage Parameter Range 


n 

2-6 

n = 2 (deg) 

0-90 

n = 3 (deg) 

0-60 

n = 4 (deg) 

0-45 

n = 5 (deg) 

0-36 

n = 6 (deg) 

0-30 

a (cm) 

1.27-5.08 

d x (cm) 

0-10.2 


The damage parameter space is sampled to obtain damage configurations, each expressed as 
a combination of input parameters (n, 0, a, cfc). The space of variables can be sampled using a 
number of different sampling methods, including random, stratified, and Latin Hypercube [40]. 
Latin Hypercube Sampling (LHS) is a type of stratified sampling method that guarantees each 
partition, or stratum, of input variable space is sampled, though not necessarily uniformly. In this 
work, LHS is performed five times for each of the variables (0, a, and d x ). Each of the five LHS runs 
corresponds to a different value of n (two, ..., six cracks) and requires the number of partitions to 
be specified. The MATLAB® implementation for LHS is used here [41], where output is provided 
in the range from zero to one. Each sample value is then scaled to the respective parameter range 
according to Table 1. 

Table 2 shows all damage configurations (26 in total) that are modeled in the ISP NN example, 
where each configuration is expressed in terms of sampled input parameters. For each damage 
configuration, the x and y dimensions of the square, damage-containing, solid FE region are provided 
in the sixth column. The x and y dimensions are each 25% larger than the diameter of the star- 
shaped damage (i.e. 1.25 * 2a), as suggested by the supplementary shell-solid boundary effect study 
described above. The last column specifies whether or not the solid, damaged region severs the 
middle stiffener. If so, the portion of the stiffener that intersects the solid model region is removed; 


10 



otherwise the stiffener remains intact. 


Table 2 Damage configurations modeled in the ISP NN example. Each damage configuration 
is assigned an alphanumeric identification with number corresponding to n. The sixth column 
provides x and y dimensions of the square region in the shell-solid ISP. 


Damage 

configuration 

a d x 6 

ID (cm) (cm) (deg) 

n Solid region Severs 

x, y dimensions (cm) stiffener? 

2A 

3.8 

7.1 

21.8 

2 

9.39 

NO 

2B 

4.2 

2.1 

87.8 

2 

10.48 

YES 

2C 

3.0 

3.6 

2.2 

2 

7.53 

YES 

2D 

1.4 

0.2 

35.7 

2 

3.49 

YES 

2E 

2.3 

9.9 

36.5 

2 

5.66 

NO 

2F 

4.5 

6.1 

44.9 

2 

11.25 

NO 

3A 

3.7 

8.3 

5.0 

3 

9.36 

NO 

3B 

1.5 

0.6 

25.2 

3 

3.72 

YES 

3C 

2.9 

9.6 

23.2 

3 

7.15 

NO 

3D 

4.0 

3.7 

32.9 

3 

9.88 

YES 

3E 

4.6 

1.7 

10.4 

3 

11.56 

YES 

3F 

2.0 

6.3 

38.6 

3 

4.92 

NO 

4A 

1.7 

6.5 

6.6 

4. 

4.15 

NO 

4B 

3.4 

1.7 

18.7 

4 

8.60 

YES 

4C 

4.9 

9.7 

25.1 

4 

12.3 

NO 

4D 

2.1 

4.4 

27.0 

4 

5.37 

NO 

4E 

3.0 

7.0 

14.9 

4 

7.52 

NO 

5A 

3.3 

8.2 

4.8 

5 

8.30 

NO 

5B 

2.2 

0.8 

6.9 

5 

5.43 

YES 

5C 

1.5 

5.4 

19.8 

5 

3.84 

NO 

5D 

3.2 

9.4 

22.6 

5 

7.91 

NO 

6A 

1.8 

8.3 

5.4 

6 

4.50 

NO 

6B 

3.7 

9.7 

26.5 

6 

9.21 

NO 

6C 

4.9 

6.0 

12.2 

6 

12.20 

NO 

6D 

4.2 

2.0 

18.4 

6 

10.61 

YES 

6E 

3.0 

3.5 

21.9 

6 

7.54 

YES 


B. Target Output: Residual Strength from Numerical Fracture Simulations 

For each input damage configuration, a numerical fracture simulation is employed to determine 
residual strength, which provides target output used to train and test the NN. FRANC3D\NG [42] 
is used to insert each parameterized star-shaped crack configuration into the solid FE region of each 
panel. An ABAQUS© contact algorithm is employed to prevent crack surfaces from overlapping 
during the applied loading. Contact properties are defined as frictionless in the tangential direction 
with “hard” pressure-overclosure behavior normal to the contacting crack surfaces, which minimizes 
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interpenetration. The FE models are then analyzed using ABAQUS®, and FE analysis results are 
post-processed to determine residual strength. 

For the sake of illustrating the NN methodology, two simplifying assumptions are made here 
to predict residual strength of the ISPs. First, the ISPs remain linear-elastic and can be analyzed 
using LEFM parameters (SIFs). Second, the residual strength can be predicted for a static crack 
configuration (i.e. crack growth is not modeled in this example). 

In the ISP NN example, the LEFM approximation of residual strength is based on mixed-mode 
I/II fracture criteria [16, 18] to account for local mode mixity (in-plane) of angled cracks in the 
star-shaped damage array. In [43], Broek described a practical mixed-mode I/II failure envelope, 
approximated by the equation of an ellipse: 


(Ki/Kj c ) 2 + ( Ku/Kuc ) 2 = 1 . ( 1 ) 

For the AA 2XXX series material in the ISP example, Kj c = 32 MPay^, and Kjj c is assumed to 
be 10% less than Ki c after results from the strain energy density criterion presented by Sih [18]. 

Using this mixed-mode LEFM-based approximation, residual strength is defined here as the 
applied traction load, Fig. 2, that first causes unstable crack growth for any point along any crack 
front of the star-shaped damage configuration. In other words, as soon as one point along one crack 
front reaches a critical combination (AT/, AT//) C on the elliptical failure envelope, the entire panel is 
assumed to fail. The method for determining the residual strength for each damaged panel is shown 
in Fig. 4 and proceeds as follows: (1) analyze the ISP FE model with P; (2) compute Kj and Kn 
at each node along each crack front using FRANC3D\NG; (3) for each crack front node , find the 
intersection point (AT/, AT//) C of the elliptical failure envelope with a straight line from the origin to 
the computed (AT/, Kjj) and subsequently find the linear scaling factor, A, that maps (AT/, AT//) to 
(Ki,Kji) c ; (4) of all the computed scaling factors, select the most critical, A c ; (5) calculate P m ax 
as P scaled by A c . 

To ensure that nonlinearity due to crack face contact does not invalidate the linear load scaling 
approach described above, each of the damaged ISPs is reanalyzed with the respective scaled load, 
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crack face contact 



Fig. 4 LEFM-based procedure for approximating residual strength of the damaged ISPs in 
the NN example. 


i.e. the approximated residual strength. In all cases, (Kj,Ku) = ( Ki,Kn) c at the predicted crack 
front failure point, indicating that the scaled loads indeed correspond to failure loads according to 
the LEFM-based failure criterion assumed for this example problem. Values of Pmax provide the 
target outputs used to train the NN. 


C. Neural Network Construction 

The inputs (sampled damage parameters) and target outputs (residual strength predictions from 
numerical fracture simulations) are used to train and test a NN. For the ISP example, a feedforward 
NN with a backpropagation learning rule [25, 44], which is a commonly used type of supervised 
NN, is constructed using MATLAB® [41]. The NN consists of a single hidden layer mapping the 
four-parameter input vectors (n, #, a, d x ) to the single- valued outputs (P ma x)- The reader is referred 
to [44] for a general discussion on NNs and details regarding specific implementation of the transfer 
functions and training algorithm described next. A tan-sigmoid transfer function is employed to 
map the weighted inputs plus bias to the interval (-1,1). A linear transfer function proportionally 
maps the weighted output plus bias from the hidden layer to the output layer. Data presented 
to the NN is divided into three sets — training, validation, and test. Weights and biases of the 
NN are adjusted at each iteration, or epoch, using the training set and a Levenberg-Marquardt 
optimization algorithm, as described in [45]. The algorithm seeks to improve performance of the 
NN by minimizing error between the NN outputs and the target outputs. Weights and biases from 
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training at any epoch are then used to check performance of the NN using the validation and test 
sets. The validation set prevents overtraining of the NN by ceasing training if performance degrades 
over a certain number of successive epochs. The test set is not used for training but is used to test 
NN accuracy following the current training epoch. The NN performance metric used here is the 
MSE, calculated as: 

1 

M ^ W = Q«E(4 8) -^) 2 > (2) 

where the superscript (i) corresponds to the training, validation, or test set, Q is the number of 
data points in the respective set, tk is target output for the k th input, and Pk is output predicted 
by the NN for the same k th input. 

The NN can be optimized by adjusting any number of parameters, including transfer functions 
between layers, number of hidden layers, various performance metrics, and so forth. In the ISP 
example, the NN is optimized by varying the number of neurons in the hidden layer (4,5,6) and by 
increasing size of the training set from 60%, to 70%, to 80% of the available data (with the balance 
equally divided between validation and test sets). Further, the performance metrics are optimized 
by minimizing MSE for the training and testing sets and by specifying that the correlation coefficient 
between NN output and targets should be at least 0.95 over the entire data set. 

D. Parametric Sensitivity Studies 

The trained NN is then employed to conduct parametric sensitivity studies, whereby sensitivity 
of residual strength to each postulated damage parameter is gauged. The sensitivity studies are 
carried out for configuration 4E, Table 2, as it represents an average damage configuration in terms 
of n, a, and d x as compared to the other configurations. 

A sensitivity study is conducted out for each of the four damage parameters. In each study, 
three damage parameters of configuration 4E are held constant while one is varied. The variable 
parameter in each study takes values in the range of the corresponding variable on which the NN 
was trained. For example, the longest a considered in the sensitivity study is no longer than the 
longest a used to train the NN, which is a = 4.9 cm in the ISP example (damage configurations 
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4C and 6C). Further, the variable parameter takes values that are equally incremented within the 
respective range. Results from the study are presented in the following subsection. 

E. Results and Discussion from Neural Network Example 

Table 3 shows the approximated residual strengths of all damaged ISPs based on numerical 
fracture analyses and LEFM assumptions outlined in II B. The table is sorted in order of increasing 
residual strength, and the corresponding damage parameters are provided to help draw preliminary 
conclusions. One immediate observation is that panels with severed stiffeners have lower (~ 50 — 
80 MPa) residual strengths, as expected. The single exception is damage configuration 2B. Though 
it severs the stiffener, configuration 2B is less critical than all other stiffener-severing cases and some 
intact-stiffener cases because it is an n = 2 configuration (straight crack) aligned with the loading 
direction. Overall, the correlation between severed stiffener and reduced residual strength highlights 
the effect of the load carrying stiffener on crack criticality. 

The ISPs with the lowest computed residual strength (configuration 3D) and highest computed 
residual strength (configuration 5C) are presented in Fig. 5. Each ISP is depicted with its respective 
residual strength, or failure load, applied. The predicted point of first-failure lies along the crack 
front indicated. Configuration 5C has more cracks and is 62.5% smaller than configuration 3D, 
though it is not the smallest of all configurations. More importantly, configuration 5C leaves the 
stiffener intact while configuration 3D results in a severed stiffener. For configuration 3D, the crack 
front that lies within the severed region and near the geometric discontinuity of the stiffener junction 
is subjected to higher stresses and is predicted to be critical. 

The optimal NN consists of four neurons in a single hidden layer with 80% of available data (i.e. 
twenty damage configurations) allocated to training. NN performance metric (MSE) as a function 
of training epochs is plotted in Fig. 6 for training, validation, and test sets. The NN is best trained 
at epoch 141, beyond which the MSE in the validation set continually increases and overtraining 
is said to occur. At this epoch, MSE of the three sets are MSE train = 0.001, MSE val = 0.87, 
and MSE test = 0.30. Weights and biases connecting the input layer (damage parameters) to the 
hidden layer and the hidden layer to the output layer (residual strength) at training epoch 141 are 
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Table 3 LEFM-based residual strength approximations for all damage configurations consid- 
ered in the ISP example, sorted by increasing residual strength. 


Damage 

configuration 

a d x Severs P max 

ID (cm) (cm) stiffener? (MPa) 

3D 

4.0 

3.7 

YES 

54.2 

6D 

4.2 

2.0 

YES 

56.8 

3E 

4.6 

1.7 

YES 

56.8 

4B 

3.4 

1.7 

YES 

58.1 

6E 

3.0 

3.5 

YES 

63.5 

2C 

3.0 

3.6 

YES 

67.6 

5B 

2.2 

0.8 

YES 

71.2 

3B 

1.5 

0.6 

YES 

73.8 

2D 

1.4 

0.2 

YES 

80.7 

6C 

4.9 

6.0 

NO 

82.7 

4C 

4.9 

9.7 

NO 

83.9 

2A 

3.8 

7.1 

NO 

89.1 

3A 

3.7 

8.3 

NO 

90.3 

2B 

4.2 

2.1 

YES 

93.8 

5A 

3.3 

8.2 

NO 

94.3 

4E 

3.0 

7.0 

NO 

98.4 

5D 

3.2 

9.4 

NO 

101.5 

2F 

4.5 

6.1 

NO 

102.1 

6B 

3.7 

9.7 

NO 

109.2 

3C 

2.9 

9.6 

NO 

110.2 

4A 

1.7 

6.5 

NO 

124.9 

3F 

2.0 

6.3 

NO 

128.2 

2E 

2.3 

9.9 

NO 

129.5 

6A 

1.8 

8.3 

NO 

130.7 

4D 

2.1 

4.4 

NO 

132.8 

5C 

1.5 

5.4 

NO 

146.7 


presented in Tables 4 and 5. 

Considering the entire set of damage configurations, the NN predictions correlate well with the 
target outputs at epoch 141, as depicted in Fig. 6. Despite the good overall correlation and small 
MSE for the training set, the MSE in the validation and testing sets (which include only three 
damage configurations each) may be too large for actual implementation onboard an aircraft. It is 
suspected that adding more samples to the entire set of damaged ISPs would further reduce these 
errors in the NN. 

The influence of each damage parameter on predicted residual strength can be visualized graph- 
ically by plotting predicted residual strength as a function of each damage parameter (see Fig. 7). 
Sensitivity can be quantified by a number of different metrics, many of which yield comparable 
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Fig. 5 Two different damaged panels (ID 5C and ID 3D) shown with respective Pmax applied. 
Panels represent damage configurations that are least critical (a) and most critical (b) of 
all configurations considered. Predicted failure point lies along the indicated crack front. 
Deformation is scaled by factor of 10. FE mesh is not shown for better contour visualization. 


Table 4 NN weights and biases used to map input layer to hidden layer for the optimized NN 
at training epoch 141. 



Input 

parameter 

Hidden layer neuron 
12 3 4 


n 

-0.43 -2.98 -1.11 1.83 

Weights 

a 

0.69 -9.17 1.61 -2.27 


X 

3.58 -0.15 1.50 -8.74 


e 

-1.53 2.87 -2.96 5.58 

Biases 


-0.42 0.45 1.37 2.94 


results [46]. Here, sensitivity is quantified by the c v , expressed as a percentage: 


c v — 


where a = 


Ttl (Prr 


- Pn 


N - 1 


(3) 


For any given sensitivity subset, i corresponds to the i th sample configuration, Pmax is the average 


17 







Fig. 6 (a) NN performance as a function of training epochs for the optimized NN; overtraining 
occurs after epoch 141. (b) Correlation between predicted and target residual strength values 
considering all damage configurations. 


Table 5 NN weights and bias used to map hidden layer to output for the optimized NN at 
training epoch 141. 



Hidden layer neuron 

Output 


1 

0.78 

Weights 

2 

0.25 


3 

-1.70 


4 

1.05 

Bias 


1.29 


residual strength of the subset, and N is the total number of damage configurations in the respective 
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subset. Sensitivities to each damage parameter are calculated as (ffl = 24.8%, c = 16.6%, 
ci n) = 6.0%, and c {0) = 1.2%. 

Orientation and number of cracks are found to have relatively minor influences on predicted 
residual strength, which is apparent both by their sensitivity metrics and by the plots (b) and (d) 
of Fig. 7. Crack length, on the other hand, has a more significant influence and causes a reduction 
in predicted residual strength as crack length increases, which is expected. What is unexpected, 
however, is the step-like behavior depicted in Fig. 7(c). This behavior is caused by binary modeling 
of the stiffener (explicitly modeling as severed or intact), a feature that is inherently implicit in 
both crack size and location. The stiffener effect is also apparent in Fig. 7(a) of damage location 
sensitivity. Predicted residual strength is lowest (and relatively insensitive to damage location) if 
the damage is located such that it severs the stiffener. As the damage location moves away from the 
stiffener and is no longer severing it, there is a linear increase in residual strength until the damage 
is located within the middle quarter of the bay. 

In general, the importance of this kind of sensitivity study is (1) to gain a better intuition of how 
and why certain damage characteristics influence residual strength and (2) to potentially decrease 
the dimensionality of the NN by neglecting parameters deemed insignificant. 

III. 3D Elastic-plastic Fracture Simulations for Improved Neural Network Training 

For discrete-source damage cases involving significant ductile tearing, a generally applicable 
3D EPFM framework should be used for improving residual strength training data. An elastic- 
plastic crack growth simulation procedure, as implemented in this section, is illustrated in Fig. 8 
and proceeds as follows: 

1. define an uncracked FE model and boundary conditions; 

2. extract a sub-region of the mesh for crack insertion, remeshing, and reconnection with the 
global mesh; 

3. map previous deformation and material state onto the remeshed model; 

4. perform nonlinear FE analysis and monitor the crack growth criterion; 
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Fig. 7 Sensitivity of predicted residual strength to each damage parameter. 


5. once criterion is satisfied, stop the current FE analysis to update crack configuration, remesh 
sub-region, and reconnect sub-region mesh with global mesh; 

6. repeat from step 3 until critical crack length is achieved or until residual strength is attained. 

The simulation procedure allows for prediction of curvilinear crack paths and arbitrary crack 
front evolution. The EPFM framework was overviewed in [47] and is described here for completeness. 
Additional details of the framework, including scripts used for implementation, are found in [48]. 


A. Nonlinear Fracture Parameter: Crack-tip Displacement 

In elastic-plastic tearing simulations, especially of thin metallic structures, crack growth should 
be characterized by an appropriate nonlinear parameter. One such parameter arises from correlation 
between crack growth and a critical amount of opening or displacement behind the crack tip (see 
[49] for details). A criterion based on this parameter, which is called the crack-tip displacement 
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I 
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ABAQUS®+Python® script 

*Map material state 
to maintain history 
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Fig. 8 Elastic-plastic crack growth simulation algorithm using FRANC3D\NG. Contributions 
from this work include evaluation of crack-tip displacement (CTD) criterion during nonlin- 
ear FE analysis and implementation of material state mapping algorithm following adaptive 
remeshing. 


[CTD) or sometimes referred to as the generalized crack (tip) opening displacement [22, 50], is 
implemented here. Notably in simulation, once a critical value, CTD cr i t , has been determined for 
a specific material and thickness through a calibration procedure, the same CTD crit is applicable 
over a range of structural configurations comprising the same material and thickness under similar 
loading. In this work, CTD is computed as: 

CTD = \jcTDj 2 + CTD// 2 + CTD,,, 2 (4a) 


CTD i = v\ — v 2 


(4b) 


CTDjj — u\ — u 2 


(4c) 


CTD n i = wi — w 2 , 


(4d) 


where u, v , and w correspond to displacements in the x, y, and 2 ) directions, respectively, and 


subscripts 1 and 2 denote the two points used to compute CTD. CTD is computed between two 


points that are initially coincident (one on each crack face) on the undeformed crack surface at 


a distance, d, behind the crack front node (i.e. in the direction normal to the crack front at the 
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particular crack front node). This is illustrated schematically in 2D in Fig. 9. In 3D, CTD values 
are computed behind multiple crack front nodes. The pair of initially coincident points where CTD 
is calculated is called a CTD point. Element shape functions are used to interpolate displacements 
(u,v,w) such that the CTD points need not correspond to nodal locations. 

y 

T* 

z 

— previous crack surface 

new crack surface 

o initially coincident points where CTD is computed 
Fig. 9 Simplified schematic of CTD implementation illustrated on crack profile. 

Crack growth occurs when CTD attains a critical value, CTD crit , within a specified tolerance. 
There are several ways to evaluate the CTD criterion when modeling a 3D crack front, including 
evaluation at a single CTD point either midway along the crack front or on the specimen’s free sur- 
face. Alternatively, the CTD criterion may be evaluated by comparing CTD crit to an average CTD 
value calculated for multiple CTD points. Because CTD crit is known to depend on 3D constraint 
at any point along a crack front, using a single CTD cr i t to predict the advance of an entire crack 
front might not be valid for all cases. A more rigorous and computationally expensive evaluation 
technique would be to compare CTD at each CTD point to a constraint-dependent CTD crit . While 
some work has been done to resolve a relationship between 3D constraint and CTD crit [20, 51], a 
constraint-dependent fracture criterion is not evaluated in the simulation described here. 

B. Material State Mapping Algorithm 

Following crack growth and remeshing, state variables are mapped from the previous mesh to 
the current mesh using an inverse isoparametric mapping routine, as in [52]. Lim et al. described 
the inverse isoparametric mapping technique for 2D elastic-plastic fracture simulations [53]. Im- 
plementation of the mapping algorithm consists of two high-level steps: (1) in the (n) mesh, state 
variables stored at integration points are extrapolated to nodes using element shape functions and 
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(2) displacements and state variables are transferred to either nodes or integration points in the 
(n + 1) mesh. The second step involves finding, for each point in the (n + 1) mesh, the natural 
coordinates (£, rj) of that point with respect to the element from the undeformed (n) mesh in which 
that point would spatially reside. The inverse problem becomes finding the natural coordinates 
(£, rf) that satisfy the known global coordinates: 

* (n+1) = 'ZN i {Z,ri)X\ n \ (5) 

where the subscript i ranges from one to the number of element nodes, X ^ are global nodal 
coordinates in the (n) mesh, X^ n+1 ) are point or nodal coordinates in the (n + 1) mesh, and N are 
element shape functions evaluated at (£, 77 ). Once (£, 77 ) are known, nodal displacements and state 
variables, C7, can be transferred from the (n) mesh to the (n + 1) mesh in a forward manner, again 
using the element shape functions, N: 


[/(«+!) =S7V i (^, n) ul n) . (6) 

Two levels of mapping are incorporated into the extended FRANC3D\NG and ABAQUS® 
software framework. First, displacements are mapped onto the undeformed mesh following crack 
growth and remeshing. Second, a mapping function available in ABAQUS® is invoked to map the 
remaining state variables (e.g. stress, strain, plastic strain) onto the deformed mesh. When mapping 
material state between successive cracked configurations, it is critical that mesh refinement in regions 
of high gradients (e.g. near crack fronts) is sufficient to minimize solution diffusion, which occurs as 
a result of extrapolation, interpolation, and nodal averaging (if employed). This effect can become 
compounded as the crack growth simulation continues. After growing the crack and remeshing, 
the updated mesh model contains additional surface area due to crack extension, and equilibrium 
must be re-established before additional load is applied. During the equilibration procedure, global 
boundaries are held fixed and the new, traction-free crack surfaces are allowed to displace in response 
to surrounding fields, as shown in Fig. 10. 
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Fig. 10 Qualitative comparison of deformation and equivalent plastic strain field after mapping 
and subsequent equilibration. Images show face of 3D mesh. Deformation is not scaled. 


C. Validation Simulation 

The EPFM framework for predicting crack growth and residual strength is validated by simu- 
lating a stable tearing test of an ISP machined from a lower wing-skin aluminum alloy, C433-T39. 
The test was conducted at Alcoa Technical Center. To illustrate the necessity of using an EPFM 
framework for predicting crack propagation and residual strength in relatively thin metallic struc- 
tures with discrete-source damage, the test is also simulated using an LEFM-based methodology. 
In the LEFM simulation, the material is modeled as linear-elastic, and crack growth is assumed to 
occur when an average value of Kj along a crack front approximately equals fracture toughness of 
the material for any increment of crack length. 

Test details, data, and results were overviewed in [54] and have also been provided to the authors 
by Alcoa Technical Center. Relevant details are described here for completeness, and additional 
details from the test program can be found in [48]. Dimensions of the panel are shown in Fig. 11(a). 
An initial two-bay saw cut of length « 2.54 cm was made at mid-height to completely sever the 
middle stiffener. The initial cut was then propagated under fatigue loading until both crack fronts 
were 2.54 cm short of reaching the intact stiffeners (2a ~ 24.1 cm). The panel was then loaded 
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monotonically in uniaxial tension until failure occurred by unstable crack growth. Crack front 
branching was observed, where an initial crack propagating toward an intact stiffener eventually 
split into two distinct cracks, one continuing into the adjacent bay and one propagating in the z 
direction within the stiffener. Photographs of the test panel with views of crack branching are 
provided in Fig. 12. 

A 3D FE model of the entire panel is constructed using ABAQUS® [1]. The FE model contains 
an initial crack of total length 24.1 cm, which corresponds to the fatigue crack length just prior to 
conducting the residual strength test. The FE model with initial crack is shown in Fig. 11(b). The 
mesh region that remains unchanged throughout the tearing simulation is modeled using 56 (731115 
and 9,400 C3D20R elements. A 38.5 x 12.7 cm 2 sub-region centered in the panel is subject to 
geometry and mesh updating within FRANC3D\NG. Depending on crack length, the sub-region 
comprises between 27,000 and 95,000 quadratic elements, including a bulk of C3D10 elements and 
a standard rosette of C3D15, C3D20, and pyramid (collapsed C3D20) elements surrounding the 
crack front (see [42] for details). The mesh interface between the sub- and global regions is coherent, 
obviating the enforcement of a coupling constraint. 

The thickened grip ends of the panel are modeled as linear-elastic with an elastic modulus 
approximately five times greater than that of C433-T39. The rest of the panel is assigned C433-T39 
material properties: E = 71.4 GPa, v = 0.3, and a y = 455 MPa [55]. The strain hardening curve 
used for C433-T39 is provided in Fig. 13. A von Mises yield criterion with isotropic hardening is 
assumed. For the LEFM simulation, the panel is modeled as linear-elastic with Kj c = 50 MPay^m 
for C433-T39 [55]. 

Boundary conditions are applied to simulate actual loading in the panel. Nodes on the bottom 
face of the lower grip end are fixed in the y direction. Displacement is applied in the y direction at 
nodes on the top face of the upper grip end. Additionally, nodes along the same top and bottom 
grip end faces are fixed in the x and z directions. For the EPFM simulation, after each increment 
of, inclusively, crack growth, remeshing, and material state mapping (see subsection III B) , all 
nodes with applied boundary conditions are held fixed while the model is brought into equilibrium 
before applying additional displacement. Additionally, based on preliminary simulation results, the 
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a) 



Fig. 11 (a) Schematic (not to scale) from [54] showing dimensions of symmetric ISP tested at 
Alcoa Technical Center. Isoparametric view of full panel indicates final fatigue crack length 
2 Oi, and cross-section view in plane of the crack shows where stiffener is completely severed. 
All dimensions are in cm. (b) Corresponding 3D FE model of ISP. Initial crack severs middle 
stiffener. Traction, P, is applied uniaxially in the y direction. 


entire back (z m i n ) face is artificially fixed after mapping and during the equilibration procedure. 
This is because resonance in the z direction is observed with increased crack growth otherwise. 
The resonance occurs when mapped tensile and compressive stresses in the faces of the panel are 
artificially reversed during equilibration of the mapped solution. The additional boundary condition 
is, however, removed after the equilibration procedure so that z displacement is allowed during the 
subsequent loading step. 

Crack growth occurs in the LEFM simulation when the average Kj value along either crack 
front reaches Kj c . A mixed-mode failure criterion is unnecessary, as Kjj and Km are negligible 
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Outboard view 

Fig. 12 Photographs of ISP with central two-bay crack from the Alcoa test program [54]. Full 
panel in load frame (left) and angled views of crack front branching into stiffener (top right) 
then exiting the stiffener (bottom right). 



Fig. 13 Strain hardening curve determined from uniaxial tension tests for C433-T39 in LT 
orientation [54]. 


(i.e. <2.5% of Ki for all crack growth increments). 

For the EPFM simulation, CTD crit was calibrated at NASA Langley Research Center from a 
middle-crack tension (MT) test of the same material (C433-T39) and thickness as the ISP [56]. In 
that work, 3D FE simulations of the MT test revealed that simulated load versus crack extension 
matched experimental data when the mode I opening angle midway along the crack front at d = 
1.02 mm reached a critical value of 6.5°. This angle corresponds to CTD crit through the relation 

tan( 6.5°) = C ^" ■ (7) 

The same criterion is applied in the EPFM simulation by specifying for both crack fronts that 
CTD cr it must attain a value of 0.116 mm at d m 1.02 mm behind the crack front and that the 
criterion be evaluated midway along either crack front. 
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Crack growth occurs in increments of ^1.15 mm (about 15% of the skin thickness), which is 
selected to be approximately the same as that implemented by Seshadri et al. in a similar simulation 
[56]. Straight crack fronts are enforced during crack growth in the skin of ISP. Upon entering the 
intact stiffener, a crack front is evolved such that (1) a realistic, arbitrary crack front profile is 
represented (though the actual evolving crack front profile was not monitored during experiment) 
and (2) the new crack front profile has relatively smooth curvature to facilitate remeshing. A cross- 
section view of the panel in Fig. 14 shows different stages of simulated crack front evolution, from 
lead crack growth in the skin, to transition crack growth within the stiffener, to complete branching. 
The simulation proceeds as depicted in Fig. 8 until both initial crack fronts completely branch and 
Pmax is attained. 



a) 


b) 



e) > f) 

Fig. 14 Cross-sectional views of ISP mesh model taken at the crack plane and magnified at 
one stiffener. A thickened red line is overlaid along the crack front (s) at each step of crack 
growth. Views show lead bay crack before entering stiffener (a); transition crack evolution 
within stiffener (b,c,d); and complete branching into two distinct crack fronts (e,f). 


Evaluation of the CTD criterion becomes nontrivial as a crack front transitions within the stiff- 
ener (i.e. while a crack front is within the stiffener but has not completely branched). Constraint 
effects introduced by the stiffener on the unsymmetric crack front profile, along with slight z dis- 
placement near the cracked region, lead to nonuniform and unsymmetric CTD values along the 
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crack front. Evaluating the CTD criterion at only one point along the crack front becomes ambigu- 
ous to implement numerically and less representative, physically, of 3D crack growth behavior. A 
simple and efficient approach to address these issues is to compare CTD cr i t to an average of CTD 
values along the crack front. For transition crack growth in the stiffener, the middle third section of 
CTD points along the crack front are averaged and evaluated to predict crack propagation. Once 
the crack fully branches, the CTD criterion is again evaluated midway along each crack front. 

D. Results and Discussion from Elastic-plastic Validation Simulation 

The consequence of using LEFM versus EPFM simulation to determine residual strength is 
made clear in Figure 15, which shows load versus crack extension for both LEFM and EPFM 
simulations. Experimental load versus crack extension was not recorded during the tests; however, 
maximum applied load is plotted for two different ISP tests of the same material and loading 
conditions. Load required to initiate crack extension is similar using either the EPFM or LEFM 
method since there are no residual stresses in the model at da= 0. Following initiation, however, the 
EPFM simulation predicts that the applied load must be increased to maintain crack propagation. 
The necessary increase in applied load occurs since a significant amount of energy in the system is 
dissipated through plastic deformation. This effect cannot be predicted by the LEFM simulation 
since plasticity effects are not modeled. As a result, much of the energy in the ISP for the LEFM 
simulation must be dissipated through creation of new fracture surface area, which means less load 
is required to drive crack growth in the LEFM simulation than in the EPFM simulation. Using the 
EPFM framework, residual strength is determined within 2% of experimental average of the two 
tests. On the other hand, the LEFM method underpredicts the average residual strength by 64%. 
Although the LEFM simulation does predict an increase in applied load at the stiffener junction 
due to geometrical effects, the increase is negligible compared to that due to plastic deformation. 

From the EPFM simulation, equivalent plastic strain field evolution in the ISP is depicted 
in Fig. 16 for the first and final crack steps. Accumulation of plastic strain in the wakes of the 
advancing crack fronts is relatively significant, extending from initial to final crack front locations. 
The general shape of 45° contour lobes at da = 0 mm is maintained at da = 44 mm both for the lead 
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Fig. 15 Applied load (traction P, Fig. 11(b), integrated over applied area) versus half-crack 
extension, da, from ISP simulation. Maximum applied load is indicated for two corresponding 
tests conducted at Alcoa Technical Center. Shaded region indicates initially intact stiffener. 


crack front extending into adjacent bay and for the crack front propagating in the z direction within 
the stiffener. At da = 44 mm, the equivalent plastic strain in each stiffener extends in the direction 
of both contour lobes to the stiffener boundary. The consistent contour lobe shapes indicate that, 
despite increased z displacement as the crack propagates and severs initially-intact stiffeners, both 
lead and branched crack fronts remain locally mode I dominant throughout tearing. 

Finally, as evident in Fig. 16 for da = 44 mm, the mapping procedure inevitably leads to imper- 
fections in the fields due to diffusion of the FE solution, which occurs in regions of high gradients 
from repeated extrapolation and interpolation procedures, see IIIB. If mapping errors significantly 
affect crack growth predictions, mesh refinement should mitigate this effect. 


IV. Conclusions 

A surrogate model methodology and 3D elastic-plastic fracture simulation toolset have been 
presented, which enable accurate residual strength prediction of damaged structures in real time. 
The methodology and toolset are particularly useful for scenarios involving metallic aircraft struc- 
tures subject to discrete-source damage during flight. An accurate prediction of structural residual 
strength in these scenarios could aid in avoidance of catastrophic crack growth and subsequent 
structural failure. 
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Fig. 16 Magnified views of simulated crack growth in ISP at half-crack extensions da = 0 mm 
(top) and da = 44 mm (bottom). Contours show evolution of equivalent plastic strain fields 
with crack growth. Deformation is not scaled. FE mesh is not shown for better contour 
visualization. Complete simulation can be viewed at www.cfg.cornell.edu . 


The surrogate model methodology relies on offline numerical fracture simulations to obtain a 
set of data points describing residual strength as a function of discrete-source damage parameters. 
Strictly for illustration, a NN has been constructed as a surrogate model for predicting residual 
strength of a representative wing sub-structure subject to discrete-source damage. In the illustration, 
offline residual strength values have been determined using computationally efficient linear-elastic 
fracture mechanics (LEFM) approximations. We have subsequently shown the consequences of 
using LEFM approximations for determining residual strength of damaged metallic structures and 
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have described an elastic-plastic fracture mechanics (EPFM) framework to accurately determine 
residual strength using high-fidelity, 3D, elastic-plastic tearing simulations. For an aluminum- alloy, 
integrally-stiffened panel exhibiting crack branching, residual strength is predicted within 2% of 
experiment using an EPFM simulation and is underpredicted by 64% using an LEFM simulation. 

The more general and rigorous elastic-plastic tearing framework should be used to generate ac- 
curate residual strength training data, especially for cases involving discrete-source damage. Also, 
the FE model for the structure of interest should include enough detail to fully capture the rela- 
tionship between a particular global loading state and onset of unstable crack growth. Furthermore, 
damage should be parameterized by taking into account onboard sensor characterization capability 
and resolution. With these considerations in mind, the general surrogate model methodology cou- 
pled with the EPFM simulation framework presented in this work provides a means of achieving 
more resilient and adaptive aircraft control. 
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